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' Foreword 


This report presents the current status and results achieved during 
the first six months of effort on Real time Flight Simulation Methodology. 
The work is a continuation of research performed during the first year 
on substitutional methods for digitization. Input signal dependent inti- 
grator approximations and digital autopilot design. Topics reported on 
herein include Digital Autopilot Design, Digital Simulation Methods for 
Linear Systems, Interactive Simulator Design Package for The Design of 
Real Time Simulators and Study of Charge Coupled Devices for Simulation. 
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I . DIGITAL AUTOPILOT DESIGN 

I .0 | nfroductlon 

In the final technical report 01 for Grant Number NASA NSG 1151, 
page 36, it was concluded: 

n The performance of a continuous autopilot which Is Imple- 
mented digitally is affected so much by the zero order hold that 
the discretization method is a secondary consideration. Tustin Is 
a relatively simple method that is satisfactory. Computer speed 
(sampling rate) may be estab I ished on the basis of the phase shift 
a designer will allow to be Introduced by the zero order hold. 

Phase lag decreases as sample speed Increases. Gain constants 
should be maintained to keep steady state errors constant. The 
relative stability of the aircraft will be decreased by the digital 
implementation so the value of computer sampling time should be 
based on the decrease in phase margin a designer is willing to 
allow. This value can be determined through sensitivity studies.” 

An example sensitivity study is presented here to demonstrate how a 
digital autopilot designer could make a decision on minimum sampling rate 
for computer specification. It consists of comparing the simulated step 
response of an existing analog autopilot and its associated aircraft dy- 
namics to the digital version operating at various sampling frequencies 
and specifying a sampling frequency that results in an acceptable change 
in relative stability. In general, the zero order hold Introduces phase 
lag which will increase overshoot and settling time. It should be noted 
that this solution is for substituting a digital autopilot for a contin- 
uous autopilot. A complete redesign could result in results which more 
closely resemble the continuous results or which conform better to original 
design goals. 

! . i Sensitivity Study 

The autopilot is an open-loop system, insofar as the complete air- 
craft system is concerned; but comparison between the analog and the digi- 
tal autopilot should be based on their performance in the overall closed- 
loop system which includes the aircraft dynamics as noted above. 
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The simplified model of an autopifot system C23 for pitch attitude 
control shown in Fig. 1.1 was chosen for the study, because It is re- 
presentative of the kind of problem a designer would typically confront.^ 
The transfer function of the short-period approximation is used to re- 
present the aircraft dynamics. An equivalent block diagram is shown in 
Fig, 1.2 In which the autopilot portion is separated from the aircraft 
transfer formation. Figure 1.3 depicts the digital autopilot system 
which Is to replace the analog autopilot of Fig. 1.2. The constants K(z) 
and H(z) of Fig. 1.3 are to be determined by applying a discretization 
method to the continuous transfer functions of the analog autopilot of 
Fig, 1.2. It Is noted that, in this example, the continuous transfer func 
tions consist of gain constants and a differentiator. By applying Euler’s 
method L33 to the differentiator*, the following discrete transfer func- 
tions for K(z) and H(z) are obtained where S(amp) = 5.6 and S(rg) = 1,19 
and K(z) Is set equal to S(amp). 


4 * 

x The transfer function representing the aircraft dynamics is for a four- 
engine jet transport flying in straight and level flight at 40,000 feet 
with a velocity of 600 ft/sec (355 knots) with the compressibility effects 
neglected. 


^Discrete differentiation operator with Euler’s method 

y(t) = u’ (t) 


UCs) 


u(t) 


Y(s ) 
y(t) 


Applying Euler’s method, yields 
u1d+UT 3 = uEnTD + Ty [nT] 


In z- domain zU(z) = U(z) + TY(z) 


Y ( z) r_-i _ Z— I 

UfzT - z dis L J T 
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H(z) = Z 


rS ( rg ) 
dts L S(amp) 


s + a = z 


r l • 19 . 

Lt-t- s + 


d Is 5.6 


□ 


= 0.212 + 


z - 


+ I 


(LI) 


The zero-order-hold and the aircraft dynamics are combined to be re' 
presented as In Fig. 1.4 where 

G(z) = Z[G, (s) D(s)3 
ho 

* Z Zl - e ~ T$ * I3.9(s t 0.306) -j 

3 s(s + I0)(s 2 + 0.805s + 1.325) 


= 13.9(1 - z"l)Zt 5 * °’ 306 a 

s 2 (s + 10) (s 2 + 0.805s + 1.325) 

= 15.0(1 ^"Mz C 0, 02509454 I °’ 059l5i363 0.0010395 f (-0.0581 54(s t 0.4025) 

s 2 s s + 10 Cs + 0. 4025 )2+l. 078422 2 


+ (-0.052705) (1.078422) y 
( s + 0.4025) 2 + I.78422 2 


= GjCz) + G 2 Cz) t G 3 (z) t Gi*(z) + G 3 (z) 


( 1 . 2 ) 


and 

G:(z) « 13.9(1 - 2-1)2^09434-, = 0,32101 I326T 

s 2 z - I 

G 2 (z) = 13.9(1 - z' 1 )Z D- , ' 059 - 1 | -' 565 ] = 0.821925946 
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G 3 (z> = 13.9(1 - z l >Zp2i20!°M^ - 


-0.01 4442 ( z - I) 
z - e"i° T 


G 4 (z) = 13.9(1 - z-ljzC-^2581 54(8+^4025 _ 

7 (s + Q.4025) 2 + I.078422 2 


( 0.808541 )(z - 1)(z - e -0.4p25T cos ] .Q7B422T) 
z 2 - 2ze"°' 4 °25 T cosl .078422T + e" 0 * 8051 


G 5 (z) = 13.9(1 - z 1 ) ZC ( ~° * 052705 ) " — j 

(s + 0.4025) 2 + 1.078422 2 

(-0.7326) (z - I) e“°‘ 42 5T sjn | ([>O7 0422T) _ 

— — ~ T ■ 0 " \ I » J) ) 

z 2 - 2ze“ 0 - 4 °25 T cosl .078422T + e "°* 805T 

in order to minimize round-off and truncation errors in the calcula- 
tions for the system response, G(z) ts put In parallel form as .shown in 
Fig. 1.5. By transforming the transfer functions of Fig. 1.5 into time- 
domain recursive relationships, the output time response can be obtained. 
Preliminary calculations using Hewlett-Packard time-sharing Basic had 
32-bit floating-point arithmetic, and the output responses obtained at 
different sampling frequencies are shown in Fig. 1.6, it Indicates that 
the word length of the Hewlett-Packard minicomputer is not large enough 
to make the truncation and the round-off errors in the calculation neglt- 
bly small. The indication appears particularly evident in Fig. 1.6 in the 
case of T = I msec (or f = 1000 cps which is the fastest sampling fre- 
quency used). 

In order to reduce calculation errors a CDC 6400 was used which has 
a word length of 60 bits, thus having 120-bit floating-point arithmetic 
in Basic. The results in Fig. 1.7 and in Table 1.1 show that the percent 
overshoot increases uniformly as the sampling frequency decreases. 












Normalized Step Responses 



Figure 1.7 Simulated Digital Autopilot System Output 
Responses Using GDC 6400. 
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II. DIGITAL SIMULATION METHODS FOR LINEAR SYSTEMS 

2.0 I ntroduction 

The methods described In this section are applicable primarily to 
linear time- invariant systems. Time varying systems whose parameters 
vary very slowly could also be treated by these techniques, and In some 
cases weakly nonlinear systems might also yield to this type treatment. 
However, the approaches are based on the solution of time invariant linear 
differential equations and, therefore., are most accurate for this class 
of problems. 

2. 1 Use of Padg Approx! mants to the Matrix Exponential for Computer 

Solutions of State Equations 

Work in this area was continued at a very low level. Previous efforts 
ClU had resulted in Fig. 2.1. This figure can be used to determine the 
percent error for a given step size and order of approx imant. Earlier ef- 

forts were concerned with negative real eigenvalues only. 

Since our last report, consideration has been given to systems with 
complex eigenvalues; however, still restricted to the left-half plane. 

It has now been determi ned that, for eigenvalues of a given magnitude, 
the percent error for the various approximants is greatest when the Ima- 
ginary part of the eigenvalue is zero. Thus, the results presented In 
Figure 2.1 can be viewed as a worst case for complex eigenvalues. The 
user need only compute the magnitudes of the various eigenvalues of the 
A matrix and define X as the largest of these. This value of X can then 
be used in conjunction with Figure 2.1. The user is guaranteed that the 
percent error actually incurred In the approximation will be no greater 
than that read off the figure. 

Tighter bounds couid be obtained by taking into account the angle of 
the eigenvalues; however, multiple plots would be required, and the re- 
sults would not be nearly as convenient for the user. 

A paper on this topic [X] has been submitted for publication to the 
IEEE Transactions on Automatic Control, 


Next, the phase margin of the digital autopilot system with different 
sampling periods was calculated and compared with that of the continuous 
autopilot system. In this calculation for which a hand calculator was 
used, up to nine digits was retained in each calculation step in order to 
ensure a degenerate accuracy. The results in Table 1.2 and Fig. 1.8 
show that, as the sampling frequency increases, the phase margin of the 
digital autopilot system converges to that of the continuous autopilot 
system as it should, 

1.2 Conclusions 

it was found that increasing the sampling period made the system 
less stable, i.e. the percent overshoot becomes larger. The phase mar- 
gin becomes smaller as expected. A designer can observe the step res- 
ponse and decide how much increased overshoot can be a! lowed. Based on 
this decision, it is simple to determine the corresponding phase margin 
from Fig. 1.9 and Its associated sampling period. 


TABLE l.l 



} 



t 

t 

i 
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Percent Overshoot vs. Sampling 

Sampling Period (sec) 

continuous case 
T = 0.0001 
T - 0.001 
T = 0.01 
T = 0.05 
1 = 0.1 


Period (Obtained with CDC 6400) 

Percent Overshoot 

27,3$ 

77.3% 

27.3$ 

28.2$ 

33.2$ 

41,0$ 


TABLE 1.2 

Phase Margin vs. Sampling Period 


Sampling Period (sec) 


continuous case 

27.4° 

T = 0.001 

27.4° 

T - 0.01 

26.7° 

T = 0.05 

24.3° 

T = 0. i 

21.0° 
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2.2 Analytical Integration of State Equations Using Interpolated Input 

2.2.0 I ntroductlon 

The basics of this approach were described in last year’s final re- 
port D3. The main idea is to use sampled values of the input signal, 
interpolate between the samples via some type of hold device and then 
use the analytical solution of the state equations to obtain the output 
signal. An important feature of this method is that the system is modeled 
exactly. The approximation is in sampling and interpolating the input. 

The solution, when one uses a first-order polynomial to represent 
the input between samples, is obtained by using the general solution 


X N ** e X N-I + / 


NT 


N- 1 


e A(NT-t) B(J(+)d+ 


(N-I)T 


and approximating U(t) by 


U( t > : U N _, ♦ (U N - v /- ( y- 'V 

The result Is 

X^ : e AT X N _! + E-A" 1 + j A“ z (e AT - l)]BU N t EA _1 C AT - i A" 2 (e AT -D^BU^ 


Regardless of the type polynomial used for the interpolation, one can 
write its transfer function, and it is simply a hold device of one form or 
another. Figure 2.2 shows a block diagram of the process. 

One would like to know the relationships between bandwidth of the in- 
put signal, W., bandwidth of the plant, W H , and the sampling frequency, W g 
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as they relate to the accuracy of the discrete representation. Certain 
minimal requirements on samp I ing frequency may be obtained by assuming 
that the hold device is perfect (that is, it has unity gain and no phase 
shift) over the frequency range 0 - “ 3/2 bu ^ does not have perfect 
rejection outside this band, and therefore higher frequencies may be 
transmitted through. Once these minimal requirements have been estab- 
lished, one can then examine the Imperfections of the hold devices over 
the frequency range 0 - “ s /2 anc * ma ^ e comparisons between the various 
ordered polynomial fits. 

2.2.1 System Bandwidth Greater Than Input Signal Bandwidth 

Consider the case where 


> to | 


( 2 . 1 ) 


u> >2w. 
s I 


( 2 . 2 ) 


Figure 2.3a shows the spectrum of the Input signal, and Fig, 2.3b shows 
the spectrum of the sampled input signal. In Fig. 3.2c we illustrate the 
characteristics of a somewhat idealized hold device with unity gain over 
the frequency range 0 - £ 1)^2 and with non-zero gain at the higher fre- 
quencies. Figure 2.3d shows the spectrum of the signal at the output of 
the hold device. In fig. 2.3e is shown an idealized gain characteristic 
of the plant. The question we wish to answer Is what must be the relation- 
ship among to , w.* and a> u so that the final output signal will be the 
same whether or not sampling occurred. From Fig. 2.3d and 2.3e it is 
clear that we must have 


0> H < (0 S - 0), 


(2.3) 


03 > (l), . , 03, 

s H + I 


(2.4) 


Note that, since 


a> H > 11) | 


(2.5) 
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the requirement of Inequality (2.4) Is more stringent than that of in- 
equality (2.2) which would be the only requirement if we assumed ideal 
low-pass fi ltering with cut-off frequency near o>| . 

Consider now the case where 


Wj_| > (flj 


and 


( 2 . 6 ) 


(d s < 2ui| (2.7) 

Figures 2.4a through 2.4e depict this situation. Looking at Figs. 2.4d 
and 2,4e, one sees that the output signals will be the same with the sampl- 
ing as it would without the sampling if, and only if, 

“H < “s " “l (2.8) 

or 


> u> H + (2.9) 

However, this contradicts inequalities (2.6) and (2.7). Therefore, we 

conclude that one cannot reproduce the output properly if 

oij_| > (2. 10) 

and 

to < 2u. (2.11) 

s I 

Thus, for the case where 


<D|_I > Uj (2.12) 

the conditions for good output reproduction are 

w s > 2tii| (2.13) 


and 


uj > tiiu + (u , (2.14) 

S n I 

Since inequality (2.14) |s the stronger, its fulfillment guarantees that 
inequality (2.13) will be fulfilled. 
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2.2.2 System Bandwidth Less Than Input Signal Bandwidth 
Now, consider the case 

W|_l < u j (2.15) 

and 

> 2oj j (2.16) 

Looking again at Figs. 2.3d and 2.3e, we see that the output will be 
reproduced properly for this case. 


Finally, consider the case 

I0|_j < U)j 


and 


(2.17) 


a) < 2a), 
s I 


Looking at Figs. 2.4d and 2.4e, we see that the output may sti 
produced properly if 


(2.18) 
be re- 


W H < “5 ~ ®l 


(2.19) 


or 


w > to,, + to. (2.20) 

S n I 

This is a less stringent condition than inequality (2.16), but it sti i i 

guarantees that the output wl II be the same with sampling as it would with- 
out the sampling. 

The final conclusion is that the same minimal conditions on sampling 
frequency 

tii w u (2.21 ) 

Sin 

apply whether 

(ii u > tii * (2.22) 

H I 

or 


(2,23) 
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These are minimal conditions and guarantee no aliasing or high fre- 
quency distortion for the case of the somewhat Idealized hold and the Ideal 
cut-off characteristics of the pi ant which we have assumed. The next step 
Is to examine the actual low frequency (0 - characteristics of the 

holds represented by the various polynomial Interpolation schemes and 
compare the accuracy taking these imperfections Into account. 

2,3 Discrete Representation Obtained via the Frequency Domain 

The design procedure developed using this approach was discussed, 
along with an application In last year’s final report D3. Since then, 
a paper K! has been written and accepted for publication in the IEEE 
Transactions on Industrial Electronics and Control Instrumentation. 
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III. INTERACTIVE SIMULATOR DESIGN PACKAGE FOR THE DESIGN 
OF REAL-TIME SIMULATORS 

3.0 I ntroductlon 

This section describes the status of research on an interactive soft- 
ware support system which will aid the design of optimum simulation models, 
The generic type of system under study is shown In Fig. 3.1. When pro- 
gramming Is completed, it is envisioned that the Simulator Design Pack- 
age (SDP) can be used to evaluate a number of different standard Integra- 
tor models (for example, Tustin, Optimum Discrete Approximation) on the 
basis of selectable error criteria or design an entirely new model suit- 
able for a particular problem. In the latter case the model would be 
designed on an interactive basis using selectable algorithms to find an 
optimal form. 

In previous work Cl] we examined a number of different substitution 
methods to determine which was most suitable under vairous error criteria. 
Based on these results, a number of substitution formulas have been chosen 
for Inclusion in SDP. Consequently, most of the work during the past six 
months has concentrated on the evaluation of optimization algorithms and 
discrete representations. At this stage, very little programming has been 
accomplished, but the elements of SDP have been defined, and programming 
is now underway. Thus, the emphasis in this report lies In Mode 2 oper- 
ation of SDP; that is, that mode of operation in which the user interacts 
with the system as it attempts to converge iteratively to an optimum model 

This section contains two sub-sections which report on different 
facets of the development of SDP. In the first sub-section we are con- 
cerned with assuming a discrete representation for a given continuous 
transfer function and then iterating to solve for optimal values of the 
parameters concerned. In the second sub-section a similar effort is re- 
ported on in which a form Is assumed for an integration operator, and 
then a random search method is used to determine the optimum parameters. 

3 . 1 Digital Simulation and Optimization via Gradient Techniques 

In previous reports we have examined the IBM, Tustin, Sage, etc,, 
methods for digital simulation. They all share a common characteristic; 
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System Configuration for Mode! Development 



once the parameters for the digital system from the continuous system are 
obtained, they cannot be altered to adjust to changes In input. These 
methods, along v/ith z-transforms, can work satisfactorily If the input 
can be categorized. However, in many physical situations, the Inputs 
cannot be categorized in this fashion; therefore, some method In which 
the digital system coefficients can be adjusted for a certain input will 
indeed be a more accurate representation of a continuous system. In 
other words, a simulation model for a given continuous system may not be 
unique but may be a function of the input. 

This model can be achieved by a variety of constraints on the res- 
ponse of the digital systems The mean square error In time domain should 
be minimum; the error in frequency domain is minimum, etc, Therefore, 
we are confronted with an optimization procedure; and the digital system 
derived by this procedure should be optimal in some way. Since most per- 
formance criteria are nonlinear, an iterative algorithm is required. This 
algorithm Is central to the process and must be modified so that Inter- 
active actions between the user and the computer are possible at al! points 
during the convergence of the algorithm. 

In order to achieve an effective optimization technique for a digi- 
tal simulation system, we must first define a performance criterion to be 
optimized, the possible constraints on the variables, the form of the di- 
gital system (so that the performance criterion can. be obtained as a func- 
tion of system parameters) and, finally, an optimization algorithm. 

3.1.1 Performance Criteria 

The natural choice here is the simulation system error, since our 
ultimate goal is to match the response of the continuous system and that 
of the digital system, i.e., the response {y^} produced by an input se- 
quence {x n > is exactly matched to the {y(nT)> sequence obtained from sampl- 
ing the continuous response. This match should occur at all points {x n K 
This exact matching is, of course, impossible; therefore, we must try to 
come as close as possible to the ideal situation. This leads to the 
definition of various errors. 
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Many error criteria have been proposed in literature. The two most 
common are the time domain error and the frequency domain error. Frequency 
domain error Is much easier to define than the time domain error, since 
we only have to compare the frequency responses of the two systems in 
the range O-T/2. Moreover, frequency responses of most systems are 
smooth; therefore, only a few comparisons are needed over the frequency 
range to obtain good results. 

During the first few months of this research, efforts have been corr* 
centrated on the frequency domain error. We will repeat the same* proce- 
dure here for 1;ie time domain error in the hope that we can obtain two 
optimization programs, one based on frequency domain error and the other 
based on time domain error. The user can try both methods for a given 
continuous system and choose the better one of the two. 

The definition of error is defined as 


E = 


M 

l 

m = I 


H(J» > - H( j« ) 
m m 


2 

, 0 <_ (til <0)2 * * * < < y 


(3.1) 


where H(jw. ) is the desired characteristic at w = <o., and H(jtu. ) is the 
simulation characteristic at the same frequency. Unlike time domain error 
definition, the frequency domain error poses some problems. First, 
by Itself, this error measure presents no constraints on poles locations. 
Using this error definition, a simulation digital system may have Its poles 
outside the unit circle In z-plane. Second, even if the poles are inside 
the unit circle, the transient responses of the simulated system and the 
digital system may be grossly different, resulting in an Inaccurate simu- 
lation. Therefore, along with this error criterion, some constraints must 
be placed on the pole locations of the digital transfer function to limit 
time domain error, as we will see later. 

3,1.2 Forms of Digital System 

The form of a digital system should be either a cascade or parallel 
combination of first-or second-order systems, because they are least 


28 


sensitive to quantization and round-off errors. Steiglitz C6] pro- 
posed the following forrru 


H(z) 


K 

l a k z 

k=l K 


- ( k- 1 ) 


N 

l 

k=l 


1 + l b k z “ k 


(3.2) . 


This form has some difficulties: first. If control is to be exercised 

over the pole locations, the denominator must be factored at certain stages 

in the minimization process; second, the pole locations may be extremely 

sensitive functions of the coefficients b . Therefore, another form Is 

K 

proposed : 


K I + a z -1 + b. z“ 2 

HU) = A n — - ( 3 . 3 ) 

k=l I + c^z" 1 + d k Z_2 

In this form there are 4k + I unknowns (A,a^, b^ ,c^,d^) to adjust for mini- 
mum error. The optimization procedure used by Steiglitz is based on the 
algorithm developed by Fletcher and PowelOH. In order to reduce the 
dimensionality, the error is minimized with respect to A, leaving 4k un- 
knowns. By supplying various starting points, the least of the minimum 
values of the error can be obtained, and, also the coefficients a^b^, 
c^,d^. Since the procedure does not place any constraints on pole loca- 
tions, after the coefficients for H(z) are obtained, we must convert ajl 
poles outside the unit circle into the unit circle and start the process 
again to obtain the final digital form. The conversion of poles (or 
zeroes) has no effect on the error. Suppose H(z) has a real pole at z 

= a. Replacing z = — is equivalent to multiplying by 

ct 



i« 


(3.4) 
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when a is on the unit circle. The same properties hold for complex poles. 
Therefore, with an adjustment In A, the definition of error Is . unaffected 
by this conversion. 


A severe disadvantage of this method Is that transient response ts 
not being taken into consideration. A pole may fall near, or on, the unit 
circle which may prevent the response from decaying as fast as that of 
the continuous response. 


Another form which was proposed by D. Schroeder DU] retains the cas- 
cade characteristic and also allows the pole magnitudes to be constrained. 
The form Is given by 


H(z) = A 


k ;, " * a 2k- I + a 2k 2 ' 2)<l + a N 2 2 ‘‘ ) 

NCP/2 ‘ NP . 

II (I + 2b 9 .cosb„. .z' 1 + b 2 z~ 2 ) n (l + b i, z “ l ) 
k=l /K ZK ~ 1 k = NCP+I 


(3.5) 


where 


NZ = number of zeroes 

NP = number of poles 

NCP = number of complex poles 

NZ 

K= largest i nteger 

(even number of zeroes) 


ef NZ - 0 if NZ = 2K 
a* NZ ^ ° i f NZ = 2K + i 


(odd number of zeroes) 


This form will be discussed in more detail later, as it proves to be the 
best form available under the frequency domain error criterin. 

3.1.3 Optimization Algorithm 

The method of Fletcher and Powell C7U is chosen here for the follow' 
Ing reasons: 

(i) it has been widely used in many similar problems reported in 
literature with a high degree of success. 


(2) It always converges, and it converges rapidly (it converges 
quadrat! ca I ly for quadratic functions). 

A number of other methods can be used, and the Fletcher-Powel I algorithm 
Is, by no means, the best ( or final) choice. Its availability as a sys- 
tem supported subroutine in some systems makes it an attractive choice 
for the Initial investigation. In the future other optimization algorithms 
will be studied In order to evaluate their performances so that the best 
possible method will be used In the proposed software package. A brief 
outline of the F I etcher-Powe II algorithm is presented here. 

The algorithm performs a one-dimensional minimization at each cycle, 
along a direction determined by the gradient and an updated estimate of 
the Hessian matrix. The generalized Taylor expansion: 

f(x + u) = f(x) + g(x)u + ~ u + G(x)u + . . . (3 6) 

where x Is an n-di mens ion vector, and u is an n-di mens Ion Incremental 
vector, g(x) is the gradient vector, and G(x) is the n x n matrix of 
second-order partial derivatives. The displacement between the point x 
and the minimum ^ is given by: 

(x mtn - x) = -G -1 (x , )g(x) (3.7) 

min mi n ® 

A method of successive linear search is used to find G -1 .' Fletcher and 
Powell suggested that initially G° be chosen as the identify matrix I; 
from this, a sequence of matrices G^ is generated which will converge 
to G” 1 * The search is terminated when the function value f(x^) has not 
decreased in the last step and the gradient vector is small. 

By itself, this algorithm Imposes no constrai nts on the variable. 
Therefore, we must place artificial limits on the pole locations before 
applying the algorithm. 

Let us now return to the minimization of error to obtain the "optima I” 
coefficients of a digital model. Recall: 
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and 


K 


H(z) = A 
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NP .,1 

JI (1+ 2b ol cos b_. -1 z" 1 + b 9l z _1 ) n (I + b.z ) 
k=l Zk Zk k = NCP+I k 


(3.9) 


The constraints on the coefficients are simple: the poles, both real and 

complex, are required to lie within a circle of predetermined radius. This 
predetermined radius is obtained from the pole locations in the s-plane. 
This restriction wifi ensure that the transient responses of the two sys- 
tems will decay at the same rate 


z-pole = e T(s ~P ,ane pole) 


(3.10) 


A less restrictive form is: 


i , „ | sT, T ( Rea I s) 

z ! 1 I® 1=0 


(3.11) 


Since the F 1 etcher-Powe 1 1 algorithm requires the function value and the 


gradient vector at each iteration, the error function must be manipulated 


ti)T 


to provide these requirements. First, replace z = aj to get 
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A can be obtained by matching the steady-state gains of the two systems: 


NCP/2 


NP 


A = H(0) 


II (l+2b 7 .cos b„, . + b 2 ) n (I + b,) 
k=l 22 K ~ NCP+I K 

^(1 +« 2k ., +a 2k XI + =' z ) 


(3.13) 


With A eliminated, we now have NZ numerator coefficients and NP denominator 
coefficients by which to minimize the error. The error can be written: 

M 

E = l CCHR - HBR ) 2 + (HI - HB1 ) 2 ] (3.14) 

i m m m m 

m- 1 

where 


HR = Real < H C j tu )) 
m m 


HI = I mag (H(jw )) 
m s 0 m 


HBR = Real (H(ju )) 
m m 


HB I m = Imag (H ( jw ) ) 


The gradient vector is given by 


known constants 


contain a. ’ s, b. 's 
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when the i th element is 
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with the - given by a similar expression. In order to obtain an exp! I- 

o D , 
i 

8HBR 3HBI 8HBR 3HBI 

cit expression for the gradient vector —g — ^ m and — — 

da. da. dDj dDj 

must be obtained. The concept Is simple, but the algebraic manipulation 
is extremely tedious. We only summarize the procedure and list the final 
results here. 

First, we use the definition of HCz) in the digital form, replacing 
A with its equivalent expression obtained from matching the d.c. gains. 
Second, we separate 

H(jw £ ) = HBR^ + JHBI ^ (A - 1,2,.. .M) (3.17) 


Third, we isolate a given variable and obtain the partial derivatives 




Rea^^hlBR^ + jHBI^)] 


(3.18) 


and 


+ JHBI t )3 

For i even: 


(3,19) 


-^-CHBR^JHBI^ 

HBR^+jHBI^ C C( * +a |„j Jcostp^T-I-aj^jCosu^Tl+Jd) 
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where 


I = Ca j | s I n&jjjT - ( H*aj_| )si n2(u a T) (3.21) 

For 1 odd; 


iI-Ehbr, + JW^ 

CHBR^+j HB I ^Xcosw^T ( I +a , + , ) - 1 -a . + , cos2ui j ) + j ( a R , s i n2u> £ T- ( I *a , + , ) s i nw g T ) ] 
[ I +a ! +a j + 1 X( 1+a { cosiu ? T+a } + ( cos2a!^T)-j (a { i nu^T+a j s i n2to A T) 

(3.22) 


For NZ odd: 


sl^HBRj, + - 


CHBR £ + jHBI^U Ctcosu^T - I) - j sinio^TU 

(I t I *F cosw^T) ~ J ®{\|Z S ^ 


(3.23) 


The derivatives with respect to the denominator coefficients are: 
For i even: 



2CHBR Jl t JHBl £ l 

+ JHBl^D .j. j + 2bjCOS bj _ i + b p 


R t j| 
R f + j I 



where 


(3.24) 


R=CCcosb, 1 +b I )(l+2b.costi3 0 Tcosb. t +b^cos2to T 

1 J "|| I X* I"- JC 
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(3.26) 


R' ~ D + 2bjCosbj_| cosus^T + b| cos 


(3.27) 


I* = “C2bj cos b._| sinus^T + b| s!n2u^TH 


(3.28) 


For t odd: 


!HBR fl + JHBI 1 = 


2CHBR a + jHBip 
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|^Rl + jlij 


(3.29) 


where 
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f l + l l & I T J Sj 


(3.32) 


1= -C2b l+J .cosb { s !nw T+b^stn2u> T] 


(3.33) 


Final ly, for i = NCP + 1, NCP + 2,..., NP 


a CHBR. + JHBlJCd - cosoi.T) + Jslnw.TU 

jl-CHBR. + JHBI. 3 = * i— 

dD ! * (I + b { )C(l + b.costo £ T) - jb .sinus jT] 


(3.34) 
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We now have explicit functions for the gradient vector, and the Fietcher- 
Powell algorithm can be applied after the constraints have been taken In- 
to account. 

Computer programs have been written for this method by D» SchroederCS] 
and at the time of this report similar local programs are being developed. 
Results will be Included In the next report. 

In summary, the following data are needed: 

(1) The form of the desired continuous system to be simulated 

TTtJw) « HR(jw) + JHI(jw) for 0 _< ai <5- <3.35) 

(2) Sampling interval T 

(3) Frequencies at which the error measure Is evaluated, t.e,, m 

n ^ 

for 0 < <oio • * • <oi < = 

1 z m T 

(4) The maximum radius of poles (to be used for constraints) 

(5) Some expected minimum error 

(6) Number of zeroes, poles, and complex poles of the digital system 

(7) Initial guesses for the numerator and denominator coefficients 

Further studies on this meThod are underway. It Is envisioned that a 
similar approach can be applied for a time domain error criterion, even 
though an accurate time domain measure, suitable for minimization techni- 
ques, appears to be complicated mathematically. 

3.1.4 Future efforts 

(1) Further investigations are to be directed in the method presented 
above. 

(2) Develop similar techniques, using time domain error criteria. 

(3) Develop general programs for IBM (which seem to be the most 
difficult from the programming point of view), Tustin and Optimum Discrete 
Approximation approach. The programming difficulties for the optimum 
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discrete approximation are at least circumvented by a study carried out 
during earlier work on this grant. It appears that In the Optimum Dis- 
crete Approximation method, a system can be broken down Into a cascade 
or parallel combination of flrst-and second-order systems for which the 
discrete approximators are readily available. Therefore, as far as pro- 
gramming Is concerned, the Optimum Discrete Approximation method can be 
considered just another substitutional method. 

3,2 Digital Simulation and Optimization via Random Search Techniques 

3.2.0 Introduction 

Real-time digital simulation of physical systems has received atten- 
tion for a number of years DO, CKO* This paper discusses a technique- 
for the development of a discrete time Integration operator to be used In 
the simulation process. The Integration operator can be optimized for a 
particular system subjected to a set of specified inputs. The class of 
systems being Investigated are those which can be represented by the set 
of state equations 

x = f(x, u) (3.36) 

where x is the n x I state vector, u-Is for the r x 1 control vector, and 
f Is the set of n functions, typically nonlinear. 

3.2.1 Integration Operator 

Figure 3.2a Is a block diagram of the mathematical relations In Eqn, 

(3.36). The vectors x and u are acted upon by the functional relations 

* 

f(x, u), prodecing the vector x. This Is then integrated to produce the 
state vector x. Figure 3.2b is a block diagram of a discrete approxima- 
tion to the continuous time system. The control vector u Is assumed to 
be sampled at a uniform rate, producing the Input samples u(k). The 
equations 

x(k) = fCx(k), u(k)H (3,37) 
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are in the same form as those representing the continuous time system. 
For example, if 

xx = cos X 2 t sin X 3 + u 


then 

* 

x^Ck) = cos xi(k) + sin x 3 Ck) + u(k) 


The function F‘(z) represents the discrete integration operator. For 
the simulation to be realizable, the denominator of F(z) must be of a high- 
er power in z than the numerator This can be shown to be correct as 

follows. Let F(z) be 


which is the familiar Tustin substitution operator for continuous integra- 


tor — . 


X(z) t 2+1 T I + z _1 

*{z) 2 z - I 2 . -{ 

I - z A 


(I - z _ 1 )x(z) = + z _ 1 )x(z) 

x(z) = z - 1 x(z) + ^-{ I + 2 “l(x(z) 
x(k) = x(k - 1 ) + j |><k) + X(k - l)H 


(3.38) 


The calculation of the state at time k is seen to depend on its derivative 
at time k. However, from the state equations, the derivative at time k is 
a function of the state at time k, (3.37), This leads to an equation of 
the form 


x(k) « x(k - 1) + I{f[x(k), u(k)3 + fCx(k - I), u(k - 1)3} 


Since the functions f(x, u) are typically nonlinear, the equation cannot 
be solved by factoring x(k) out of the right-hand side of the equation. 
Therefore, (3.38) represents an unrealizable simulation. If the integra- 
tion operator had been chosen as 


F(z) 


T,z+I 
~ 2 z - I 


•) z 



then 


x(k) = x(k - I ) t ^ Cx(k - 1 ) t *(k - 2)3 (3.39) 

which is a realizable simulation. Adding the z" 1 term to the Tustin opera- 
tor, however, degrades its performance. It has been shown, though, that 
to be closed-loop realizable, the power of the denominator must exceed 
that of the numerator. What is desired is a closed- loop realizable opera- 
tor which can be optimized for a particular system being driven by a set 
of known inputs. 

To this end, a discrete-time integration opera’^r of the following 
form is chosen. 

N . 

T . l y 

F(z) - 4 ~ 0 (3.40) 


where T is the sampling period, and the X’s are a set of free parameters, 
the values of which are to be optimized. This operator yields a realizable 
simulation, since the power of the denominator is always one greater than 
that of the numerator. The pole at z = I corresponds to a pole at the 

4_h 

origin in the complex s-plane, and the N order pole at the origin corres- 
i'h 

ponds to an N 1 order pole at infinity in the s-plane Cl G 
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The free parameters tn F(z) are optimized using an Idealized model 
form of a model reference adaptive control system Cl 23. Figure 3,3 Is 
a block diagram of this form. An exact, or nearly exact, solution to 
the system’s differential equations represents the ideal model of the 
system* The discrete simulation, including the parameters to be optimized, 
represents the controlled process, A set of inputs Is applied to the 
model and the process, and the error at each sample time Is squared and 
summed. At the end of one run, the free parameters are perturbed under 
the control of an optimization technique; and the sequence Is repeated. 

This continues until the error of the digital simulation Is sufficiently 
sma 1 1 . 


3.2.2 Optimization Technique 

The perturbation of parameters in F(z) is controlled by an Adaptive 
Random Search Optimization technique (ARSO). Random perturbation methods 
have been shown to solve a large class of optimization problems faster 
than gradient techniques when the number of unknown parameters exceeds 
four D33, Cl 43, Cl 53- In addition, the convergence time has empirically 
been shown to increase linearly with the number of unknown parameters 
rather than exponentially or quadratical iy [J53* 

In the ARSO technique, the mean and variance of a uniformly distributed 
random variable are adaptively selected for each unknown parameter based 
on recent successful experiments. In this way the step size and direction 
are controlled by past successful perturbations. Adaptive step size has 
been shown by Schumer and Stelglltz Cl53 to be a powerful technique in 
multidimensional problems without ridges or valleys. Adding adaptive 
step direction should tend to broaden the range of applicability. 

The performance index used with ARSO is a vector valued one; that Is, 
one which requires simultanelous minimization of all components Cl33. 

This allows consideration of several cost functions, such as integral 
square error, minimum energy, etc., at the same time. In the nonlinear 
problem to be studied, air-craft dynamics, the mean square error for each 
of the state variables Is used as a component in the performance Index. 
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That is, 



N 

I tx 
J « I U 



and 


(3.4!) 


j = Di J 2 ••• J n 3T (3.42) 

where N is the number of sample periods per trial, x. . is the calculated 

*t*h fh ^ J 

value or the i Tn state variable at the j sample time, and x|\ Is the 

exact value. For a trial to be considered a success, at least one component 

of J must be reduced, and no component may Increase in value. 

The unknown parameters are perturbed in the following manner. 

X.(j t |) = X* t 6X . (j) (3.43) 

fh 

where X ( j + I) Is the new value of the I parameter, X* Is the ''best- to- 

JiU * 

date" value of the i T parameter; that is, the value of X. when the mlnlmum- 

to-date value of the J vector was calculated, and 6X.(j) Is the random 

i*h ^ 

perturbation for the i parameter. The perturbation is calculated as 
fo 1 1 ows : 

6X } (j) = Vj(J) + v / 3o2(jT [2 RND(o) - ll (3.44) 

f h 

where y,(j) is + he current value for the mean of the i random variable, 
a 2 (j) is the current value for the variance, and RNO(o) is a uniformly 
distributed random variab le on the interval C'0, O* Equation (3.44) pro- 
duces a random number from a uniformly distributed random variable with 
mean u and variance a 2 . At the present time, the variance - - r-.i same 
for each of the parameters, while a mean is calculated for each. The use 
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of Individual variances, as well. Is an area to be studied, Stability 
considerations may place constraints on the parameter calculated In (3.43). 
If a particular perturbation places a value outside Its limit for stabi- 
lity, the value may be moved deterministically Inside the limit. The par- 
ticular system will determine what the stability limits for the parameters 
are. 

When a particular trial Is successful; that is, 

J I <_ I £ i <_ n 

+h 

where Jy Is the minimum value of the i element of J, the means of the 
random variables are updated. 

U j (j ) = xf - nrij <3.45) 

where 

c,(l) c. (2) c. (3) c. (4) c. (5) 

m . - -V + -V- + -V- + -V- + -nr- (3 - 46) 

The c.’s are past values of A|; that is, previous values of the "best-to- 
date" parameters. c (I) is the most recent value of X|. Since the most 
recent best value corresponds to a smaller performance index than a pre- 
vious best value, (3.45) tends to select the most favorable direction for 
the next perturbation. For an unsuccessful trial, the mean is not changed 
until the number of consecutive failures becomes large. The strategy for 
the mean and variance is as fol lows: 

Initially, the mean Is zero; and the variance Is equal to two. This 
permits large perturbations In the values of the parameters. As the trials 
become successful, the means are computed by (3.45); and the variance Is 
greatly reduced. This permits the parameter values to fol low favorab le 
terrain. The mean is updated at each success as previously mentioned. 

When the number of consecutive failures reaches 100, the variance is in- 
creased to allow larger step sizes. At this point, the means remain 
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unchanged. If the number of consecutive failures reaches 1000, the 
variance Is greatly increased; and the means are set to zero* This allows 
the possibility of Jumping out of a local minimum to an area with a small- 
er performance index. When the number of consecutive fat lures reaches 
10,000, the program Is terminated. If, at any time, e successful trial 
occurs, the variance is reduced; and the means are updated as before. 

3.2.3 Preliminary Results 

Thus far, the ARS0 technique has only been used for the devopment of 
general-purpose Integration operators for linear systems. A three-element 
performance index was used to measure deviation of the Integrator operators 
from Ideal Integration. Frequency domain magnitude error from zero to 
one-half the sampling frequency, phase error over the same frequency range, 
and phase error from zero to eight-tenths of one-half the sampling fre- 
quency made up the three components of J. One and four parameters were 
used in F(z) on different experiments. The two operators obtained are 
shown below. 

F(z) 

F(z) » T(.97l8z 3 - . I768z 2 - .0386z + .1613) 

z 3 (z - I ) 


The performance index for these two operators, as well as that for the 
forward difference operator, which is closed-loop realizable, are shown 
in Table 3.1. The operators were also evaluated In the time domain by 
using them to solve a system of seven linear differential equations using 
a fourth-order Runge-Kutta solution as the Ideal model. Table 3.2 sum- 
marizes the results. The performance Index was of ’the form in (3.41) 
and (3.42). Reasonable agreement between frequency domain optimization 
and time domain, evaluation is realized for the linear system. This is 
expected from application of Parseval’s Theorem |J63* 
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TABLE 3.1 

Frequency Domain Optimization 
RMS ERRORS 


Operator 


ARSO - I 
ARSO - 4 
Forward Diff. 


1.18161 dB 
2.61323 dB 
I .76163 dB 


52.9347 degrees 
50.5897 degrees 
52.9347 degrees 


42.5422 degrees 
36.5964 degrees 
42.5422 degrees 


TABLE 3.2 

Time Domain Evaluation 
RMS ERRORS 


ARSO - 


.217174 

4.99025 

53.9662 

319.485 

1004.97 

1683.85 

1034.9 


.53871 
1 1 .4281 
109.522 
574.329 
1550.95 
1872.31 
1014.4 


Fwd. 

Diff. 


,531912 

1 1 . 

.5529 

M2, 

.056 

583, 

.429 

1498. 

.52 

1482 

.75 

661 

.88 
















/ f2 ( + )d + / F(s) p { „ s)ds (3.49) 

—CD v ™*C0 

which relates the frequency and time domains. Of course the user should 
formulate his performance index to reflect what the major concern is. 

The following transfer function was put In state variable form and used 
for the evaluation. 

x 36 S 2 + 2.31 + 2.72 s + 1.65 s 2 + 7.25s + 81 

PCs) = — — — * * • ~ * — — - 

s 2 + 8.4s +36 s 2 + 5.62s +3.! s + .62 l.!25s z + 13. 33s + 81 

(3.50) 

This is the transfer function for an autopilot and therefore is of interest 
to us. 

3.2.4 Further Study 

The optimization technique described here will be applied to simulat- 
ing the dynamic behavior of an aircraft. Initially, motion In a vertical 
plane will be considered, with more degrees of freedom added as familiarity 
with the optimization technique is gained. Modifications in the strategy 
of ARSO, such as individual variances for the random variables, will also 
be studied. 



IV. STUDY OF CHARGE COUPLED DEVICES FOR SIMULATION 
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4.D Introduction 

The objective of this work was to investigate the feasibility of us- 
ing CCDs (charge coupled devices) In real-time aircraft simulation. When 
a single digital computer is used for simulation, the number of calcula- 
tions made serially Is so large that even with relatively large computers 
the sampling rate cannot be made fast enough to make discretization error 
negligibly small. For this reason CCD devices were primarily examined 
to determine If the sampling rate could be made large. 

Charge coupled devices are relatively new. The Beil Telephone Labora- 
tories announced the initial discovery of CCDs In May cf 1970 Cl 7- 193. 

A CCD Is a semiconductor component C203 which is capable of storing quan- 
tities of charge and moving them practically intact from one storage loca- 
tion to another. The quantities of charge represent the magnitude of 
signal samples. If the movement takes place at regular time intervals, 
a discrete delay line action is achieved. The addition of signal extrac- 
tion, weighting, and other manipulations such as summing yield analog 
data processing of untque potential. Obviously, CCDs can aiso be used 
to process digital data. 

The CCDs belong to the class of components called charge transfer 
devices of which the "bucket brigade device" (BBD) C2l3 Is a member. The 
BBD is older but of less potential value for signal processing than the 
CCD C223. 

For an in-depth review there are several papers available which des- 
cribe the physical principles of CCDs, how they are constructed, and de- 
tails about their operational performance C233. Briefly, they consist of 
a substrate of doped semiconductor materia! (such as n-type) on which a 
very thin layer of insulation is placed. Cells are then formed by deposit- 
ing metal electrodes on the insulating surface. When the electrodes are 
biased with respect to the substrate, a potential well is formed in the 
semiconductor materials. Charge is injected into the device with a PN 
junction (in the case of n-type substrates). The charge may then be moved 
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between adjacent cells by appropriately pulsing the associated electrodes. 
After transfer of charge from a cel I, the electrode voltage may be re- 
turned to a quiescent value. The transfer between cells Is not perfect 
but is quite adequate for many applications, having an efficiency on 
the order of 0.9999 to .99999, 1 .e. the transfer inefficiency e is on the 
order of I0“ 4 to I0~ s [24D. Thus, a single packet of charge can be moved 
through many cell locations without obtaining significant signal distor- 
tion. The e depends on clock frequency as one would expect. At tne 
present state of development, operation is usually restricted to leas 
than [0 MHz. Transfer efficiency Is important to analog data processing 
because, unlike digital processing, signal levels may not be re-established. 
The thermal generation of eiectron-hole pairs is a serious source of dis- 
tortion, especially at elevated temperatures. This restricts storage in 
a cell te times less than one second (usual ly, one or two orders of magni- 
tude less than one second) at room temperature. Storage time is reduced 
by a factor of two for every 10°C rise in temperature. This is illustrated 
later with calculations. 

The types of operations performed include ail kinds cf linear trans- 
formations such as correlations, matched filtering, discrete Fourier and 
Hadamard transforms, Karhounen-Loeve expansions, etc. in addition, fast 
time simulation can make use of discrete convolution. With some peripheral 
circuitry, these operations may be implemented by various combinations of 
the basic building blocks C253 consisting of: 

1. Serial In/Serial Out (SI/SO) Registers 

2. Serial In/Paral lei Out (SI/PO) Registers 

3. Parallel In/Seria! Out (Pi/SO) Registers 

The SI /SO Register is the simplest form of CCD processor. It con- 
sists of a large linear array of cells forming a discrete analog shift 
register or delay line. This device is shown in Figure 4.1a. 

As shown in Figure 4.1b, the basic SI/PO structure can be used to 
extract desired data points after a given delay time. If tap weights are 
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used wi+h current summing on the outputs, then it is possible to imple- 
ment various linear transformations. For example. Figure 4.2 depicts 
a non-recursi ve or transversal filter whose output at the nth clock 
period is the weighted sum of the previous k samples, 

k 

g(n) = y h.x . , n > k (4,1) 

tV " - 1 - 

This is, of course, discrete convolution. If the tap weights h. In (4.1) 
define samples of the time inverse of a desired signal, then the CCD pro- 
cessor of Figure 4.2 is a matched filter. 

The PI/SO building block is shown in Figure 4.1c. Here, the inputs 
are synchronously sampled and stored in particular cells of the linear 
array. 

4.1 Simulation Considerations 

The most appropriate approach for the simulation being considered is 
the discrete convolution described by (4.1) which applies to linear statton' 
ary systems. Unfortunately, the equations governing aircraft systems are 
both non-stat ionary and nonlinear; but for operating conditions where the 
assumption of small perturbations is valid, the equations can be linearized 
at each equilibrium condition of interest. This removes the linearity 
and time-varying restrictions but at the price of either having variable 
tap weight devices or a large number of devices which can be switched 
in as required. Practical variable tap weight devices are not available, 
and development progress appears uncertain so the second approach was 
assumed. 

The idea here is to use a bank of CCD transversal filters, as shown 
in Figure 4.3, each of which covers one specific equilibrium condition of 
the system to simulated in order to cover the whole operating range de- 
sired in the aircraft simulation. 

For small perturbations the motion of an aircraft can be considered 
to consist of modes with relatively long periods (such as phugoid and 
spiral divergence) and modes with relatively short periods (such as shbrt- 
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period and roll subsidence). As will be discussed later, there is a 
practical limit on the total delay of CCD transversal filters, e.g., 

100 m sec at room temperature. Even though much longer delays can be 
achieved by cooling the devices total delay of more than 10 to 20 sec 
is not desirable because discretization effect becomes significant due 
to the required decrease in sampling frequency. This is due to the 
limited number of stages of a CCD transversal filter. Therefore, the 
first compromise will be to use a small computer to digitally simulate 
modes with long periods while using CCD transversal filters for modes 
with short periods. 

4.2 Error Sources 

The transfer inefficiency e per unit cell represents the fractional 

portion of the signal charge which is left behind when a transfer takes 

place. At the present, CCDs are constructed with typical values of e ~ 

IQ-k-to |Q“ 5 . The e remains approximately constant below clock frequencies 

of f =1 MHz. Above this e starts to increase due to the finite transit 
c 

time of the signal charge from one cell to another. The e also depends 
on the quantity of charge in a cell. A d-c bias change called a "fat- 
zero” is required to make e less dependent of signal amplitude. The 
transfer inefficiency affects the CCD operation such that the input signal 
is dispersed or spread out in time, thus limiting the filter length, i .e. 
the number of delay stages. 

The mechanism for dark current (or leakage current) is the thermal 
generation of electron-hole pairs. It causes the stroage wells to be 
slowly filled with minority carriers which gradually mask the stored In- 
formation. Since the amount of charge is proportional to time, the effect 
is to limit the low frequency operation of CCD devices. It is found ex- 
perimentally that total delays up to 100 m sec at room temperature can be 
achieved before dark current effects become significant C26]. This Is 
true for. a 1 1 devices, irrespective of size or number of elements. However, 
the dark current is temperature dependence and decreases by a factor of 
approximately 2 for every I0°C drop In temperature. Thus, dark current 
can be minimized by suf f icient coot I ng of the CCD. This would be prac- 
tical for the simulation of large systems. 
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Tap weight characteristics are a source of simulation error. There 
is a practical limit on the accuracy of the weighting coefficients in CCD 
transversal filters. There is also a limitation of the range of maximum 
to minimum values imposed by practical consideration. This requires that 
tap weights smaller than the minimum realizable value either be set to 
the minimum realizable value or zero. 

Because the CCD transversal filter has a finite number of elements 
and the sampling frequency is limited, the total delay that can be achieved 
with CCD transversal filters Is also limited. Thus, non-zero values in 
a weighting function at instants of time exceeding the total delay have 
to be truncated. Maximum charge holding capacity, i .e. saturation, can 
Introduce simulation errors. Errors due to linearization have not been 
studied 

4.3 Rough Estimate of Required Number of Filters 

The number of transversal filters required is a very important para- 
meter and Is considered in the following. Generally, the coefficients 
of linearized dynamic equations of motion for aircraft systems are deter- 
mined by the following ' f light conditions: 

( 1 ) Mach number 

(2) Angle of attack 

(3) Altitude 

(4) c.g. location 

Suppose the si x-degree-of-f reedom equations of motion are uncoupled. 
Then, a typical longitudinal set of equations would have four different 
excitation inputs: (a) elevator deflection, (b) flap deflection, 

(c) throttle change, and (d) spoiler deflection. Since longitudinal 
motion has three-degree-of-f reedom and each degree-of-f reedom is Imple- 
mented with a separate set of CCD devices, a total of twelve CCD trans- 
vers I filter is necessary for each different equi I ibrium condition of 
the longitudinal motion. That means that, in case of coupled motion, 
approximately eighty CCD devices would be necessary for each different 
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equilibrium condition for six-degree-freedom simulations. Considering 
that the four flight conditions can vary independently to cause changes 
in the coefficients of the equations of motion and considering the wide 
range in stability derivatives, particularly due to Mach number, altitude 
and angle of attack, it is not difficult to imagine that thousands of CCD 
transversal filters would be needed to adequately cover the whole range 
of interest in ordinary aircraft simulation, 

4.3 Future Work 

The study was not persued in depth because the first rough estimates 
of performance for modes with large settling time and the number of CCD 
section required for realistic flight operation indicated that present 
technology favors digital technology. Use of microprocessors in parallel 
operation seems far more promising now than CCDs, so the research has 
been moved in that direction. Fundamental architecture and software pro- 
blems are being considered to try to achieve real time simulation that is 
superior to present large computer systems. 

The use of CCDs for a single set of flight conditions may be practi- 
cal but this question was not addressed because microcomputers offer a 
possible solution to the more general problem. 


Appendix 4. I 


Dark Current Considerations 

As stated in the text, the achievable total delay of CCD transversal 
filters is limited, due to dark current. At room temperatue, total de- 
lay up to 100 msec can be achieved before dark current effects become 
significant In many applications. 

In order to increase the useful storage time, cooiing of CCD devices 
have been considered. In the following It is shown that, by cooling a 
CCD device from room temperature (300°K) to 240° K, the dark current is 
reduced approximately 500 times below than that of room temperature. 

That means that, theoretically, total delay as large as 50 seconds can 
be achieved at 240 D K before dark current effects become significant. 

Dark current has six possible sources C273: 


(1) Generation of carriers via generation-recombination centers 
In the depletion region. 

(2) Generation of carriers via surface states at the S i — S i O 2 . 
i nterface. 


(3) Avalanching at the channel-stop boundary. 

(4) Diffusion of minority carriers out of the neutral region of 
the bulk. 


(5) Processing errors. 

(6) Tunneling between bands. 


Among the above six possible sources of dark current, carrier genera- 
tion in the depletion region (Source I) and carrier generation via sur- 
face states at the si ! icon- si I icon-dioxide interface (Source 2) are the 
dominant sources C27]. Thus, the dark current, J Q , can be wrltlen as 


J 


D 


qn,W qn.S 
^3 ,10 

2t 2 


- n r (|^t 



(4.2) 
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where 


rij = intrinsic carrier concentration 
W = depletion region width 
t = bulk I Ifetime 

S - surface recombination velocity when the surface is depleted, 
o 

q = charge of an electron 
Us I ng the f ormu I a 


S 0 = N ss irkT V th " 


(4.3) 


where 


N = the number of surface states per cm 2 per eV 

*3 b 

T = temperature (°k) 

a = capture cross-section for electrons or holes, whichever is 
srna 1 1 er 

V_j.^ - /3kT7m = thermal velocity of the electrons (4.4) 

t = C/V„ (4.5) 

th 

with 


C = temperature independent constant 
yields 

I? “ v th ° ^ < 4 - s) 

and 

qS 

y«T’V th «T 3 / 2 (4.7) 

The ratios of these values at two different temperatures are 
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i 


X 


^1t ) T - 240°k .. ^240 _ n 0 
- 0.9 


iSSL 

1 2t't * 300°k 


(4.8) 




(--A) 


2 J T = 240 °k _ ,240, 3 ^ 2 _ „ 
c " l 300 J “ U * 


qS 

(_£.) 

2 'T = 300° k 
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(4.9) 


Thus, 


,qW qS o } 

0.72 < -il— -g-T = 240 ° k < 0.9 


( qW + qS o } 

v 2t ~ 2~*7 = 300°k 


(4.10) 


From (4.10) it is observed that the term + —A varies little with 


'2t ’ 2 


temperature. 

Now, calculate the change of n. with temperature change. It is 


n. 2 = A TV Eg/kT 
I o 


C4.il) 


where - constant 
Energy Gap 

Boitzman’s constant (=8.62 x IO _5 eV/°K) 


o 

eg 

k 


Thus, 


(n 2 \ 

i J T = 240° K _ 


(y3) { E ~^9/kT, 

T = 240® K ■ T = 240 °K 


Cn 2 ) 


r t3 \ 

i 'T = 300® K T =- -srirtO! 


'T = 300° K 


-Eg/kT. 


(4. I2> 


(e 


T = 300® K 


SO 
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(4.13) 


<T3 ’t “ ^ . ( ^°) 3 . Q . 5.2 

(T 3) 500 

11 ; T « 300°K 


and 


(Eg/kT).j. _ 240°K 


1.12 eV 


(8.62 x KT 5 eV/°K) (300°K) 


= 55. i 


( Eg/kDj _ 3oo° K 


1.12 eV 


(8.62 x I 0“ 5 eV/°K) (300°K) 


= 43.31 


and, it follows that 


(4.14) 


( e -Eg/kT )T ^ = £ -55.i = I.J76 x j o -2 ** 


(e -Eg/kT) T ^ = £-43.31 « | ,55| x !0" 19 


(4.15) 


Therefore, 


(n j 2) T = 240°K 
(n j 2) y = 300°K 


0.512 


- 3.88 x 10 

1 .551 x 10 19 


6 


(4.16) 


and 


p T« 240°K .. |>97 x l0 -3 , <5001-1 (4.17) 

' n l' T = 300°K 

From this calculation it Is found theoretically that dark current becomes 
500 times smal ler at T = 240°K than that at T = 300°K. 
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